Electron Transfer Rates in Polar and Non-Polar Environments: a Generalization of Marcus’ Theory to Include an Effective Treatment of Tunneling Effects

A multistep kinetic model in which solvent motion is treated in the framework of Marcus theory and the rates of the elementary electron transfer step are evaluated at full quantum mechanical level is proposed and applied to the calculation of the rates of intramolecular electron transfer reactions in rigidly spaced D–Br–A (D = 1,1′-biphenyl radical anion, Br = androstane) compounds, for five acceptors (A) in three organic solvents with different polarity. The calculated rates agree well with experimental ones, and their temperature dependence is almost quantitatively reproduced.

E lectron transfer (ET) reactions are ubiquitous in chemistry and play a relevant role in modern technological applications such as solar cells, molecular electronic devices, biosensing, and photocatalysis. Marcus' theory is the cornerstone for understanding mechanistic aspects of ET reactions and still represents the most robust theoretical framework for the rationale design of new chromophores to employ in technological devices. 1−3 Indeed, it relates the kinetics of the ET process to physically sound quantities: the ΔG 0 of the whole ET reaction and the reorganization energy associated with the motions of both solvent and reactants. However, because of the classical treatment of nuclear motion, Marcus' rate expression is often unable to reproduce the observed temperature dependence of ET rates, with valuable exceptions where, however, the temperature dependence of the reorganization energy and the reaction free energy is accounted for, 4 posing problems for a deeper understanding of important mechanistic aspects of several phenomena, as for instance charge transport in organic semiconductors, where tunneling likely plays an important role. Further problems have been found for strongly exothermic ET processes, for which predicted rates are many orders of magnitude slower than the observed ones. 5 There are several outstanding works which have extended Marcus' theory to include quantum effects due to high frequency modes, 6−13 but, to the best of our knowledge, the task of including the whole heat bath provided by intramolecular coordinates of the redox pair, which was revealed to be crucial for reproducing the unusual temperature dependence of ET from bacteriopheophytin anions to primary ubiquinone in bacterial reaction centers, 14−16 has not been undertaken yet for systems in which polar solvents are expected to play a significant role.
With that purpose in mind, herein, we generalize a recently proposed multistep kinetic model of ET reactions, 17,18 which by separating the motion of the solvent from the internal one of the redox pair, makes it possible the use of a treatment of tunnelling effects which includes the whole set of nuclear coordinates of the redox pair and accounts for changes of both the equilibrium nuclear positions and vibrational frequencies upon electronic transition, as well as the effects due to normal mode mixing. 14,19,20 The proposed model is then successfully applied to a set of electron transfer reactions in rigidly spaced donor-bridge-acceptor systems, for which ET rates and, in some cases, their temperature dependence have been experimentally determined in solvents of different polarity. 21,22 The multistep mechanism is shown in Scheme 1, where D − A and DA − denote the initial and the final states, respectively, and {D − A}* and {DA − }* denote the ensembles of transient structures in vibronic resonance with each other.
Step 1 is an activation step which brings the donor and the acceptor species into electronic degeneracy, thus triggering step 2, i.e., the elementary ET reaction. The backward reaction of step 1 and the forward one of step 3 account for relaxation of all of the nonequilibrium species to their minimum energy structures, including also the environment. The activation step is necessary in all of those case in which the initial and final states are not in vibronic resonance at frozen solvent coordinates, a situation which can often occur in polar environments, because the solvent polarization free energy can be larger than the ΔG 0 of the ET reaction. For the activation step, we thus assume that only solvent coordinates are involved and that its rate constant has the usual Arrhenius dependence on temperature: i where k 0 is a transmission coefficient and ΔG # is the standard activation free energy, i.e., the free energy difference of {D − A}* with respect to D − A, which, closely following Marcus' reasoning, but at frozen intramolecular coordinates, 3 is given by where ΔG fi 0 is the free energy difference between initial and final states, λ s and λ i are the reorganization energies of solvent and of the redox pair A and D, respectively, and λ = λ s + λ i . Derivation of eq 2 is straightforward and is based on the same assumptions made by Marcus in his original work, 1,3 but for the sake of clarity it is reported in the Supporting Information. It is noteworthy that ΔG # is higher than in Marcus' theory, because of our assumption that only solvent coordinates are involved in the activation step.
Backward step 1 and forward step 3 consist of the solvent response to a nonequilibrium charge distribution of the solute. Those processes have been extensively studied by time dependent spectroscopic measurements (Stokes shifts), and their rates have been experimentally determined. 23,24 Experimental values from time dependent spectroscopic measurements will be used throughout for k dact . The latter quantity is related to k 0 by the principle of microscopic reversibility, which must hold for the activation step, even though neither equilibrium nor steady state has been assumed, because of its purely mechanical origin, based on the invariance of the mechanical equations of motion under the operation of time reversal. 25,26 The principle of microscopic reversibility implies the following relation between k act and k dact , averaged over all of the possible paths bringing degenerate reactants to degenerate products: 25 where g act and g eq are the total degeneracies of the activated and the equilibrium states. The above condition implies: Since the degeneracy factors g i cannot be easily estimated, here we have left k 0 as the only solvent dependent adjustable parameters to be fixed by comparison with experimental rates.
Step 2 is the elementary ET process under resonant conditions, whose rate can be reliably calculated by resorting to first order time dependent perturbation theory, i.e., the Fermi Golden Rule (FGR), according to which the rate of a nonradiative transition between two electronic states |i⟩ and |f⟩ is given by ing integration over nuclear coordinates. 19 From f(λ), ρ(ΔE,T) can be obtained by using the inverse Laplace transformation. 14,19,20,28,32 Very accurate estimates of ρ(ΔE,T) of both radiative and nonradiative transitions have been obtained in the past by using this procedure. 35−39 The chemical structures of the donor and the acceptors considered here are drawn in Figure 1. The energies of the initial A D − A and final DA − states, computed at the density functional theory (DFT) level, using the continuum polarizable medium (PCM) approach for evaluating the contribution of solvent polarization, see Computational Details, are reported in Table 1, for the three solvents, iso-octane, tetrahydrofuran (THF), and dibutylether (DBE), for which experimental measurements are available. 21,40 Calculated electronic energy differences for ET from donor to each acceptor agree well with the corresponding experimental ET free energies (ΔG 0 ), as shown in Figure S1, where theoretical ΔE fi are plotted against measured ΔG 0 's. Hereafter, we will consider the computed ΔE fi as free energies, assuming, as usual, 3 that the entropy of reaction due to intramolecular vibrations is negligible. Entropy contributions of solvent motions are instead accounted for by λ s .
All ET reactions are exothermic when both solvent (Q) and molecular (q⃗ ) nuclear coordinates are in their equilibrium conditions, but in the present kinetic model the elementary ET step occurs at fixed Q, either at Q 0i , the equilibrium solvent coordinate of the initial state, at when ΔG eff 0 < 0, or at Q c , the point at which the potential energy surfaces of the initial and final states cross each other when intramolecular coordinates are kept fixed at their initial equilibrium value. Thus, to apply the proposed kinetic model, it is first necessary to ascertain whether or not ET requires solvent activation, which in turns requires the evaluation of λ s . The latter quantity can be estimated by using several approaches: the Born-Onsager model, 1,41 molecular dynamics simulations, 42,43 the nonequilibrium polarizable continuum model, 44 and microscopic theory of solvent response. 45,46 Herein, taking advantage of the availability of experimental ET rates at different ΔG, λ s 's are directly estimated from the experimental rate constants of the extremely exothermic BIP-BQO or BIP-NQO pairs and, from the electronic coupling elements reported in ref 43, upon assumption, which can be easily verified a posteriori, that those ET processes do not require solvent activation, i.e. ΔG eff 0 < 0. The adopted approach has the limitation that λ s does not depend on the specific   Parson,43 in those systems, λ s is weakly dependent on the acceptor. The above procedure, see the caption of Figure 2 for more details, yields λ s = 0.75 eV and λ s = 0.53 eV for THF and DBE, respectively, obtained as averages between the λ s of the two exothermic BIP-BQO and BIP-NQO pairs (0.86 and 0.63 eV in THF and 0.58 and 0.47 eV in DBE, for BIP-BQO and BIP-NQO, respectively). Estimated solvent reorganization energies are in reasonable agreement with those obtained by Parson using quantum mechanics/molecular dynamics (QMMD) simulations 43 and compare very well with those estimated by Closs and co-workers, fitting experimental data in THF with a Marcus like rate expression with a single high frequency quantum mode. 21 Estimated λ s 's are plotted against solvent dipole moments in Figure S2; a remarkable linear dependence has been found.
For the nonpolar iso-octane solvent, the above procedure yields λ s 's vanishingly small for all redox pairs, as it occurs in photosynthetic reaction centers, 47 so that in iso-octane, ΔG eff 0 = ΔG fi 0 . All of the calculated rate constants in iso-octane are reported in Table 2, together with experimental ones. For NQO, the calculated ET rate significantly disagrees with the experimental value, but, as already noted in previous papers, 21,43,48 this strongly exergonic ET process could also involve an electronically excited state. Indeed, if the ET rate is evaluated at ΔG fi 0 corresponding to the first excited states of the NQO anion, obtained by TDDFT computations, keeping fixed the electronic coupling element, a reasonable agreement with the experimental value is found, see Table 2. According to TDDFT results, also ET from BIP to BQO could involve the first electronically excited state of BQO, but in this case better agreement with the experimental rate is obtained for ET into the ground state, see Table 2.
For most of the A/D systems in THF and DBE, ΔG eff 0 is positive, see Table 1, so that ET requires an activation step. ET rates are therefore determined using Pauli's master equation approach, solving the system of differential equations associated with the kinetic scheme of Scheme 1, with k act and k ET calculated from eqs 1−5, and k dact taken from experimental time dependent spectroscopic measurements: k dact = 7 × 10 11 and k dact = 8 × 10 10 , for THF and DBE, respectively. 23,49 The k 0 of eq 1 has been taken as a solvent dependent adjustable parameter; satisfying agreement with experimental results has been obtained by setting k 0 = 5 × 10 12 and k 0 = 4 × 10 11 for THF and DBE, respectively.
Calculated ET rate constants in THF and DBE are reported in Table 2 together with experimental ones, and in Figures 3  and 4, where theoretical and experimental rate constants are reported as a function of ΔG of ET reactions. According to   22,48 For the BIP/BQO pair, ΔG eff 0 < 0, and ET does not require a solvent activation step, whereas ET in the BIP/NAP pair needs activation by solvent motion. The calculated rates at different T are reported in Figure 5, which clearly highlights the differences between the two ET processes: the temperature dependence of ET from BIP to NAP is dominated by the Arrenhius exponential factor, which makes ET rates change by over 4 orders of magnitude by changing T from 180 to 400 K. ET from BIP to BQO occurs by tunneling, and since the temperature dependence of ρ(ΔE, T) is modest, ET is almost independent of temperature. Although our treatment has neglected important factors, such as the temperature dependence of ΔG eff 0 and λ s , 50 which proved to be crucial for simulating the unusual temperature dependence of ET rates observed in the fullerene-porphyrin dyad, 4 the proposed approach provides a satisfying picture of the T dependence of ET. It is noteworthy that the dielectric constant of THF exhibits a significant temperature dependence, which, in the range −78 to 30°, can be expressed by the following relation: ϵ(T) = −1.50 + 2650/T. 51 Since ϵ(T) increases as T decreases, λ s will also increase, increasing the ET rates, see Figure 2. Thus, a more refined treatment, such as those outlined in refs 45 and 46, could further improve the agreement with experimental data.
We have presented a simple generalization of Marcus' theory, which allows the inclusion of the whole heat bath provided by intramolecular coordinates of the redox pair in the evaluation of the tunnelling contribution to ET rates, a crucial task for properly handling vibronic degeneracy. In the proposed approach, the elementary ET step occurs always by tunneling, which can be or not be triggered by the environment, according to energy conditions: for highly   The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter exothermic ET processes, for which the ΔG 0 of ET overcomes the reorganization energy of the environment, the initial and final states are in vibronic resonance, and ET can occur by tunneling at frozen environmental coordinates. Vice versa, when the energy stabilization of a polar environment is larger than the ET free energy variation, environmental motion is required to trigger ET tunneling, and the whole process is described by a multistep mechanism, in which also the rates of solvent stabilization of a nonequilibrium charge distribution play an important role. The faithful reproduction of the ET temperature dependence in two different cases testifies to the versatitility of the proposed mechanism; of course, further applications are needed for a better assessment of its general applicability.
Most of this work has been stimulated by a series of recent papers by Parson, who used QMMD simulations to compute ET rates. 43,52−54 Our approach avoids the use of QMMD simulations, relying on the separated treatment of the nuclear motion of the environment from that of the redox pair, whose interplay is treated, when necessary, at the Pauli master equation level. That is not a great limitation; cases in which coherent effects play a role have been successfully addressed using Pauli master equations in the past. 18 In principle, the approach could represent a physically sound alternative to QMMD simulations, and indeed a significant improvement has been found, if enough experimental data are available, since two quantities, the solvent reorganization energies and k 0 , have to be evaluated from experimental data; in practice, when experimental data are scarce, it is still necessary to find a valid approach to evaluate from first-principles λ s and k 0 . Work is in progress along this line. The approach also provides a significant improvement with respect to nonadiabatic theories developed in the past, for the most exothermic ET reactions in a nonpolar solvent. 21 Although the procedure is entirely numerical, the approach offers a mechanistically very clear picture of the ET process. The proposed multistep mechanism also shares common facets with other recently published mechanistic pictures: the flickering resonance 55 and the unfurling mechanisms 56 are both based on the assumption that transient degeneracy among different redox species triggers coherent charge motion.

■ COMPUTATIONAL METHODS
Franck−Condon weighted densities of states have been computed by using a development version of the MolFC package, available on request. The internal (curvilinear) representation of normal coordinates has been adopted in all of the cases. 57 Equilibrium geometries, normal coordinates, and vibrational frequencies in the neutral and anionic forms were computed at the density functional theory (DFT) level by using the B3LYP functional with the 6-31+G(d,p) basis set. Excited state energies have been obtained by time dependent DFT computations, using the same functional and basis set. In all of the computations, the solvent effect was included by using the equilibrium polarizable continuum model (PCM). The electronic coupling elements reported in ref 43

■ ACKNOWLEDGMENTS
We are indebted to Prof. W. Parson for having suggested this work and for many useful discussions. Financial support of University of Salerno is gratefully acknowledged.